In [1]:
import sys
!{sys.executable} -m pip install -r requirements.txt
Requirement already satisfied: cvxpy==1.0.3 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from -r requirements.txt (line 1))
Requirement already satisfied: numpy==1.13.3 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from -r requirements.txt (line 2))
Requirement already satisfied: pandas==0.21.1 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from -r requirements.txt (line 3))
Requirement already satisfied: plotly==2.2.3 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from -r requirements.txt (line 4))
Requirement already satisfied: scipy==1.0.0 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from -r requirements.txt (line 5))
Requirement already satisfied: osqp in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: scs>=1.1.3 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: six in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: fastcache in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: ecos>=2 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: toolz in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: multiprocess in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: pytz>=2011k in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from pandas==0.21.1->-r requirements.txt (line 3))
Requirement already satisfied: python-dateutil>=2 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from pandas==0.21.1->-r requirements.txt (line 3))
Requirement already satisfied: requests in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: nbformat>=4.2 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: decorator>=4.0.6 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: future in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from osqp->cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: dill>=0.2.6 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from multiprocess->cvxpy==1.0.3->-r requirements.txt (line 1))
Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: idna<2.7,>=2.5 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: urllib3<1.23,>=1.21.1 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: certifi>=2017.4.17 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from requests->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: ipython_genutils in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: traitlets>=4.1 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: jsonschema!=2.5.0,>=2.4 in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4))
Requirement already satisfied: jupyter_core in /Users/ehamel/anaconda3/lib/python3.6/site-packages (from nbformat>=4.2->plotly==2.2.3->-r requirements.txt (line 4))
You are using pip version 9.0.1, however version 18.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.
In [2]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import cvxpy as cvx
import plotly as py
import plotly.graph_objs as go
import helper

py.offline.init_notebook_mode(connected=True)
%matplotlib inline
plt.style.use('ggplot')
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
RuntimeError: module compiled against API version 0xc but this version of numpy is 0xb

Optimization with an Alpha Model and a Risk Model

In this exercise, we will construct an optimization problem using a stock universe consisting of 3 stocks. This problem will inherently be somewhat artificial, but using a stock universe of 3 stocks will allow us to create plots in 3 dimensions. These are useful for illustrating exactly what is going on in each step and will help develop your intuition for the problem. In this exercise, we will do the following steps:

  1. Create a 1-year-momentum-based alpha vector from some stock price data.
  2. Create a risk model using PCA.
  3. Construct the optimization problem using the alpha vector and risk model, and apply a set of standard constraints.
  4. Solve the optimization problem.

Load the data

Load the data from the file stock_prices_advanced_optimization.csv.

In [3]:
prices = pd.read_csv('stock_prices_advanced_optimization.csv', parse_dates=['date'], index_col=0)

The data are the price trends for 3 stocks, Stock A, Stock B and Stock C, for 4 years, from 2013 to 2017.

In [4]:
prices.head()
Out[4]:
A B C
date
2013-07-02 100.000000 102.000000 104.000000
2013-07-03 99.872139 102.298629 104.821312
2013-07-05 99.109975 102.733517 104.256297
2013-07-08 99.971278 102.877047 103.631199
2013-07-09 98.866192 102.568696 106.375325
In [5]:
prices.plot()
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x113cadc88>

Calculate the returns of these price data.

In [6]:
rets = prices.pct_change()[1:].fillna(0)

Create an alpha vector

In [7]:
from scipy.stats import zscore

We'll create an alpha vector based on 1-year momentum.

In [8]:
# 1-yr momentum alpha

def log_returns(series):
    return np.log(series[-1])-np.log(series[0])
    
alpha = prices.rolling(window=252).apply(log_returns).rank(axis='columns').apply(zscore, axis='columns')

Now we'll take the row of most recent alpha values to use as our alpha vector. There should be a value in the vector for each stock.

In [9]:
# Take most recent set of values
alpha_vector = alpha.iloc[-1,:]
print(alpha_vector)
A    0.000000
B    1.224745
C   -1.224745
Name: 2017-06-30 00:00:00, dtype: float64

For this problem, we'll use the optimization objective $-\boldsymbol{\alpha}^T \mathbf{x}$. Remember, we are trying to minimize this function (to maximize alpha). Let's plot $-\boldsymbol{\alpha}^T \mathbf{x}$ as a function of the components of $\mathbf{x}$ so that we can get a better sense of the "shape" of the function.

In [10]:
n = 10
x = y = z = np.linspace(-2,2,n)

xv, yv, zv = np.meshgrid(x, y, z, indexing='ij')
obj_val = np.full(xv.shape, np.nan)

for i in range(n):
    for j in range(n):
        for k in range(n):
            obj_val[i,j,k] = -alpha_vector[0]*xv[i,j,k]-alpha_vector[1]*yv[i,j,k]-alpha_vector[2]*zv[i,j,k]
In [11]:
hover_text = helper.generate_hover_text(xv, yv, zv, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C', fcn_values=obj_val, fcn_label='Objective Function')

trace1 = go.Scatter3d(
    x=xv.flatten(),
    y=yv.flatten(),
    z=zv.flatten(),
    text = hover_text.flatten(),
    mode='markers',
    marker=dict(
        size=4,
        color=obj_val.flatten(),     # set color to an array/list of desired values
        colorscale='Viridis',   # choose a colorscale
        colorbar=dict(
                title='Objective Function'
            ),
        opacity=0.4,
        showscale=True
    ),
    hoverinfo = 'text'
)

data = [trace1]
layout = helper.create_standard_layout('Alpha Vector Optimization Objective Function', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

In this plot, the value of the objective function is represented by the color at each point in weight space. Notice that the function increases with weight on Stock C, decreases with weight on Stock B, and doesn't depend on the weight on Stock A.

What happens if we try to run the optimization now, by trying to minimize this objective function, under no constraints?

In [12]:
def find_optimal_holdings(
    alpha_vector):
    
    # define weight vector
    x = cvx.Variable(len(alpha_vector))
    
    obj = cvx.Minimize(
                -alpha_vector.values.flatten()*x
            )
    
    constraints = [] # no constraints
    
    prob = cvx.Problem(obj, constraints)
    prob.solve(verbose=True, max_iters=500)

    optimal_weights = np.asarray(x.value).flatten()

    return pd.DataFrame(data=optimal_weights)
    
In [13]:
optimal_weights = find_optimal_holdings(
    alpha_vector
)
print("Optimal weights: ", optimal_weights)
ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -0.000e+00  +0e+00  0e+00  1e+00  1e+00  1e+00    ---    ---    0  0  - |  -  - 
 1  -8.918e+01  -0.000e+00  +0e+00  0e+00  1e+00  9e+01  1e-02  0.9890  1e-04   0  0  0 |  0  0

UNBOUNDED (within feastol=0.0e+00).
Runtime: 0.000743 seconds.

Optimal weights:        0
0  None

As expected, without constraints, the problem is unbounded. We could achieve a smaller and smaller objective function value by taking larger and larger long positions on Stock C, and larger and larger short positions on Stock B. But this would increase risk and leverage. This is where our constraints come into play.

Create a Risk Model

We are going to create a PCA-based risk model for our 3 stocks. First, let's plot the returns data and their mean, just to have a sense for the shape of the data.

In [14]:
m = rets.mean()
In [15]:
# Trace for mean return vector m
mean_vec = np.vstack((np.full(3, 0), m.values)).T

hover_text2 = helper.generate_hover_text(mean_vec[0], mean_vec[1], mean_vec[2], 'Mean of Returns of Stock A', 'Mean of Returns of Stock B', 'Mean of Returns of Stock C', sig_figs = 7)

trace2 = go.Scatter3d(
    x=mean_vec[0],
    y=mean_vec[1],
    z=mean_vec[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='red',
        opacity=0.9,
        symbol="diamond"
    ),
    name = 'mean daily return',
    text = hover_text2.flatten(),
    hoverinfo = 'text'
)


# Trace for data
hover_text3 = helper.generate_hover_text(rets['A'].values, rets['B'].values, rets['C'].values, 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace3 = go.Scatter3d(
    x=rets['A'].values,
    y=rets['B'].values,
    z=rets['C'].values,
    mode='markers',
    marker=dict(
        size=4,
        color='#7FB3D5',  
        opacity=0.3,
    ),
    name = 'daily returns',
    text = hover_text3.flatten(),
    hoverinfo = 'text'
)

data = [trace2, trace3]

layout = helper.create_standard_layout('Returns Data', 'Return of')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

There are a few things to notice here.

  1. The vector of mean returns is very small compared to the volatility of the data! Admittedly, this is an artificial example with fake data, but this is also generally the case.
  2. The volatility of Stock C is larger than the volatility of Stock B, which is larger than the volatility of Stock A.

Given these data, we can already set up some expectations for what the results of PCA will look like. Remember that we expect the first PC to point in the direction of maximum variability of the data. Given the volatility of Stock C, we expect the first PC to point approximately parallel to the "Stock C axis". If we get results that are very different from this, we know that something is wrong with what we are doing. Let's now run the PCA algorithm, and see what we get.

First we'll mean-center the data.

In [16]:
rets = rets - m

Then we'll create a class that will fit our risk model.

In [17]:
from sklearn.decomposition import PCA

class RiskModelPCA():
    
    ANN_FACTOR = 252
    
    def __init__(self, num_factors):
        self._num_factors = num_factors
        self.num_stocks_ = None
        self.factor_betas_ = None
        self.factor_returns_ = None
        self.common_returns_ = None
        self.residuals_ = None
        self.factor_cov_matrix_ = None
        self.idio_var_matrix_ = None
        self.explained_variance_ratio_ = None

    def fit(self, returns):
        self.num_stocks_ = len(returns.columns)
        mod = PCA(n_components=self._num_factors, svd_solver='full')
        mod.fit(returns)
        
        self.factor_betas_ = pd.DataFrame(
            data=mod.components_.T,
            index=returns.columns
        )
        
        self.factor_returns_ = pd.DataFrame(
            data=mod.transform(rets),
            index=returns.index
        )
        
        self.explained_variance_ratio_ = mod.explained_variance_ratio_
        
        self.common_returns_ = pd.DataFrame(
            data=np.dot(self.factor_returns_, self.factor_betas_.T),
            index=returns.index
        )
        self.common_returns_.columns = returns.columns
        
        self.residuals_ = (returns - self.common_returns_)
        
        self.factor_cov_matrix_ = np.diag(
            self.factor_returns_.var(axis=0, ddof=1)*RiskModelPCA.ANN_FACTOR
        )
        
        self.idio_var_matrix_ = pd.DataFrame(
            data=np.diag(np.var(self.residuals_))*RiskModelPCA.ANN_FACTOR,
            index=returns.columns
        )
        
        self.idio_var_vector_ = pd.DataFrame(
            data=np.diag(self.idio_var_matrix_.values),
            index=returns.columns
        )
        
        self.idio_var_matrix_.columns = index=returns.columns

    def get_factor_exposures(self, weights):
        B = self.factor_betas_.loc[weights.index]
        return B.T.dot(weights)

Let's fit the risk model with 2 factors (i.e., we'll keep 2 PCs).

In [18]:
rm = RiskModelPCA(2) # create an instance of the class with 2 factors
rm.fit(rets) # fit the model on our data

Let's look at the factors (PCs).

In [19]:
rm.factor_betas_
Out[19]:
0 1
A -0.015530 0.059674
B -0.055588 -0.996725
C -0.998333 0.054570

The first PC points almost completely in the "Stock C" direction. The second points almost completely in the "Stock B" direction. This makes sense given the fact that we see most variability in the "Stock C" direction, and if we collapse that variability, we'd see the next most variability in the "Stock B" direction. Let's plot the PC vectors to make this more clear.

In [20]:
PC_scaler = 0.04 # The PC vectors have length 1, but this is larger than the scale of our data, so for visualization purposes, let's plot scaled-down versions of the PCs. 

# Trace for PC 0
pc0 = np.vstack((np.full(3, 0), rm.factor_betas_[0].values)).T*PC_scaler

hover_text4 = helper.generate_hover_text(pc0[0], pc0[1], pc0[2], 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace4 = go.Scatter3d(
    x=pc0[0],
    y=pc0[1],
    z=pc0[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='#45B39D',
        opacity=0.9,
        symbol="diamond"
    ),
    line=dict(
        color='#45B39D',
        width=3
    ),
    name = 'PC 0',
    text = hover_text4.flatten(),
    hoverinfo = 'text'

)

# Trace for PC 1
pc1 = np.vstack((np.full(3, 0), rm.factor_betas_[1].values)).T*PC_scaler

hover_text5 = helper.generate_hover_text(pc1[0], pc1[1], pc1[2], 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace5 = go.Scatter3d(
    x=pc1[0],
    y=pc1[1],
    z=pc1[2],
    mode='lines+markers',
    marker=dict(
        size=4,
        color='#FFC300',
        opacity=0.9,
        symbol="diamond"
    ),
    line=dict(
        color='#FFC300',
        width=3
    ),
    name = 'PC 1',
    text = hover_text5.flatten(),
    hoverinfo = 'text'

)

# Trace for data
hover_text6 = helper.generate_hover_text(rets['A'].values, rets['B'].values, rets['C'].values, 'Return of Stock A', 'Return of Stock B', 'Return of Stock C')

trace6 = go.Scatter3d(
    x=rets['A'].values,
    y=rets['B'].values,
    z=rets['C'].values,
    mode='markers',
    marker=dict(
        size=4,
        color='#7FB3D5', 
        opacity=0.3,
    ),
    name = 'daily returns',
    text = hover_text6.flatten(),
    hoverinfo = 'text'
)

data = [trace4, trace5, trace6]

layout = helper.create_standard_layout('Returns Data with Factor (PC) Directions', 'Return of')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

Let's look at the fraction of variance explained by the two factors we kept. This should confirm our intuition about the relative scale of variance explained by the first two factors.

In [21]:
plt.bar(np.arange(2), rm.explained_variance_ratio_, color=['#45B39D', '#FFC300']);
plt.title('% of Variance Explained by Each Factor');

Now let's look at the factor returns. We'll convert them to "price series" and visualize the evolution of these over time. Recall that these are the data expressed in the "factor" (PC) basis––the projections of the data vectors onto the factor dimensions.

In [22]:
(rm.factor_returns_ + 1).cumprod().plot(color=['#45B39D', '#FFC300'])
Out[22]:
<matplotlib.axes._subplots.AxesSubplot at 0x11a256860>

Note that in this example, since the first factor is so close to the direction of the return on Stock C, the factor return (converted to a "price series") for factor 0 looks a lot like the price evolution of Stock C, but inverted. This is because the PCs are vectors that represent directions in space—they are only unique up to a sign change, and thus the projections onto them are unique up to a sign change.

In [23]:
prices.plot()
Out[23]:
<matplotlib.axes._subplots.AxesSubplot at 0x11a2f3898>

Create Optimization Constraints

As we saw above, we need to constrain our optimization problem. Let's now create and plot the following constraints:

  • risk, based on our risk model
  • leverage
  • market neutral
  • factor exposure limits
  • individual weight limits

Risk

In [24]:
B = rm.factor_betas_
F = rm.factor_cov_matrix_
S = np.diag(rm.idio_var_vector_.values.flatten())

Using the $\mathbf{B}$, $\mathbf{F}$, and $\mathbf{S}$ matrices, let's write a function that takes in the portfolio weights and spits out the portfolio risk. Remember that the portfolio risk is calculated $\mathbf{x}^T(\mathbf{B}^T\mathbf{F}\mathbf{B} + \mathbf{S})\mathbf{x}$, but be careful to ensure the matrices are all oriented correctly.

In [25]:
def risk_fcn(x):
    f = np.dot(B.values.T, x)
    return np.dot(f.T, np.dot(F, f)) + np.dot(x, np.dot(S, x))
In [26]:
n_samples = 25
grid_max = 2.5
risk_grid, spacing, xv, yv, zv = helper.evaluate_fcn_on_grid(grid_max, n_samples, risk_fcn)

Now let's look at a plot of the value of the portfolio risk as a function of the portfolio weights in 3 dimensions.

In [27]:
hover_text = helper.generate_hover_text(xv, yv, zv, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C', fcn_values=risk_grid, fcn_label='Portfolio Risk')

trace7 = go.Scatter3d(
    x=xv.flatten(),
    y=yv.flatten(),
    z=zv.flatten(),
    mode='markers',
    marker=dict(
        size=4,
        color=risk_grid.flatten(), 
        colorscale='Viridis',
        opacity=0.3,
        showscale=True
    ),
    text = hover_text.flatten(),
    hoverinfo = 'text'
)

data = [trace7]
layout = helper.create_standard_layout('Portfolio Risk as Modeled by our PCA Risk Model', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

This plot uses a color at several points in 3-D space to represent the portfolio risk as a function of weight on each stock. Note that it increases with weight on each stock. This helps us visualize the risk function, but what does it look like to constrain risk? Let's visualize the boundary of the region within which risk is below a limiting value. This will help us visualize the shape of the space we search with optimization. We will do this for each of our constraints.

In [28]:
risk_data = helper.plot_iso_surface(risk_grid, 0.05, 25, 2.5, "Portfolio Risk = 0.05", '#F5B7B1', True)

Leverage

Now let's write a function that takes in the portfolio weights and spits out the portfolio leverage, $\sum_i|x_i|$. This will allow us to visualize the space that satisfies the constraint $\sum_i|x_i| \leq 1$.

In [29]:
def lev_fcn(x):
    return np.abs(x[0]) + np.abs(x[1]) + np.abs(x[2])

n_samples = 25
grid_max = 2.5

lev_grid, _, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, lev_fcn)
lev_data = helper.plot_iso_surface(lev_grid, 1.0, n_samples, grid_max, "Leverage Ratio = 1", '#1F618D', True)

Market Neutral

For the market neutral constraint, we constrain the sum of the weights. Let's write a function that calculates this quantity, given the portfolio weights. Then we can visualize the plane $\sum_i x_i = 0$.

In [30]:
def mn_fcn(x):
    return sum(x)

n_samples = 25
grid_max = 2.5

mn_grid, _, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, mn_fcn)
mn_data = helper.plot_iso_surface(mn_grid, 0, n_samples, grid_max, "Sum of Weights = 0", '#D35400', True)

Factor Exposures

Let's also calculate the factor exposures. Then we'll plot the planes defined by the factor exposure limits. The optimization will constrain the solution to lie between each pair of planes. We'll constrain each factor exposure to be between -0.1 and 0.1.

In [31]:
def fac_fcn(x):
    return np.dot(B.values.T, x)

grid_max = 2.5
n_samples = 25

fac_grid, spacing, _, _, _ = helper.evaluate_fcn_on_grid(grid_max, n_samples, fac_fcn)

factor_limit_data = []

for factor in range(B.shape[1]):
    for l in range(2):
        mult_fac = l*2-1
        factor_limit_data.extend(helper.plot_iso_surface(fac_grid[:,:,:,factor], 0.1*mult_fac, n_samples, grid_max, 'Factor Exposure Limits', '#D2B4DE', False))

layout = helper.create_standard_layout('Factor Exposure Limits', 'Weight on')

fig = go.Figure(data=factor_limit_data, layout=layout)         
        
py.offline.iplot(fig)

Individual Weight Limits

Finally, let's look at the space inside the constraints on each individual weight. We'll constrain each individual weight to be between -0.55 and 0.55.

In [32]:
x_max = 0.55

x = np.array([-1, -1, 1, 1, -1, -1, 1, 1])*x_max
y = np.array([-1, 1, 1, -1, -1, 1, 1, -1])*x_max
z = np.array([-1, -1, -1, -1, 1, 1, 1, 1])*x_max
hover_text = helper.generate_hover_text(x, y, z, 'Weight on Stock A', 'Weight on Stock B', 'Weight on Stock C')


weight_data = [
    go.Mesh3d(
        x = x,
        y = y,
        z = z,
        colorscale = '#FCF3CF',
        intensity = np.full(8, 1),
        i = [7, 0, 0, 0, 4, 4, 6, 6, 4, 0, 3, 2],
        j = [3, 4, 1, 2, 5, 6, 5, 2, 0, 1, 6, 3],
        k = [0, 7, 2, 3, 6, 7, 1, 1, 5, 5, 7, 6],
        name='Weight on Stock B',
        showscale=False,
        opacity = 0.3,
        hoverinfo='none'
    )
]

trace = go.Scatter3d(
    x=x,
    y=y,
    z=z,
    mode='markers',
    marker=dict(
        size=6,
        opacity=0.0001,
        color='#BFB1A8', 
    ),
    text = hover_text.flatten(),
    hoverinfo = 'text',
    showlegend=False
)
    
weight_data = [weight_data[0], trace]    
    
layout = helper.create_standard_layout('Individual Weight Limits', 'Weight on')
fig = go.Figure(data=weight_data, layout=layout)
py.offline.iplot(fig)

Plot all the constraints on the same axes

Finally, let's visualize the space that satisfies all the constraints. It has to be inside all of the enclosed spaces, on the plane defined by the market neutral constraint, and between the planes defined by the factor exposure constraints.

In [33]:
data = risk_data
data.extend(lev_data)
data.extend(mn_data)
data.extend(factor_limit_data)
data.extend(weight_data)
layout = helper.create_standard_layout('Visualize Intersection of All Constraints', 'Weight on')
fig = go.Figure(data=data, layout=layout)
py.offline.iplot(fig)

Optimization

Finally, let's define and run the optimization problem with the same objective function as above, but this time, with all the constraints we just discussed. Remember that the optimization is looking for the point in the space that satisfies all the constraints that minimizes the alpha function.

In [34]:
def find_optimal_holdings(
    alpha_vector,
    risk_model,
    risk_cap=0.05,
    factor_max=10.0,
    factor_min=-10.0,
    x_max=0.55,
    x_min=-0.55):
    
    B = risk_model.factor_betas_
    F = risk_model.factor_cov_matrix_
    S = np.diag(risk_model.idio_var_vector_.values.flatten())

    x = cvx.Variable(len(alpha_vector))
    f = B.values.T*x
    
    risk = cvx.quad_form(f, F) + cvx.quad_form(x, S)

    obj = cvx.Minimize(-alpha_vector.values.flatten()*x)

    constraints = [
        sum(cvx.abs(x)) <= 1.0,
        sum(x) == 0.0,
        x <= x_max,
        x >= x_min,
        risk <= risk_cap,
        B.values.T*x <= factor_max,
        B.values.T*x >= factor_min
    ]
    
    prob = cvx.Problem(obj, constraints)
    prob.solve(verbose=True, max_iters=500)

    optimal_weights = np.asarray(x.value).flatten()

    return pd.DataFrame(data=optimal_weights, index=alpha_vector.index)
In [35]:
optimal_weights = find_optimal_holdings(
    alpha_vector,
    rm
)
ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +4.621e-18  -5.915e+01  +1e+02  1e-01  7e-01  1e+00  3e+00    ---    ---    1  1  - |  -  - 
 1  -4.203e-01  -5.106e+00  +1e+01  9e-03  1e-01  7e-02  3e-01  0.9071  5e-03   1  1  1 |  0  0
 2  -1.079e+00  -1.892e+00  +2e+00  1e-03  2e-02  2e-02  5e-02  0.8543  4e-02   1  1  1 |  0  0
 3  -1.228e+00  -1.334e+00  +2e-01  2e-04  3e-03  2e-03  6e-03  0.8750  6e-03   1  1  1 |  0  0
 4  -1.224e+00  -1.227e+00  +6e-03  5e-06  7e-05  8e-05  2e-04  0.9811  8e-03   1  1  1 |  0  0
 5  -1.225e+00  -1.225e+00  +6e-05  6e-08  7e-07  9e-07  2e-06  0.9890  1e-04   1  1  1 |  0  0
 6  -1.225e+00  -1.225e+00  +7e-07  6e-10  8e-09  1e-08  2e-08  0.9890  1e-04   1  1  1 |  0  0
 7  -1.225e+00  -1.225e+00  +8e-09  7e-12  9e-11  1e-10  2e-10  0.9890  1e-04   1  0  0 |  0  0

OPTIMAL (within feastol=9.1e-11, reltol=6.4e-09, abstol=7.9e-09).
Runtime: 0.000539 seconds.

In [36]:
optimal_weights
Out[36]:
0
A 3.857385e-11
B 5.000000e-01
C -5.000000e-01

The optimal weights are 0 weight on Stock A, 50% long on Stock B and 50% short on Stock C.

In [ ]: